{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Regression: Using a Decomposition (Cholesky Method)\n",
    "--------------------------------\n",
    "\n",
    "This script will use TensorFlow's function, `tf.cholesky()` to decompose our design matrix and solve for the parameter matrix from linear regression.\n",
    "\n",
    "For linear regression we are given the system $A \\cdot x = y$.  Here, $A$ is our design matrix, $x$ is our parameter matrix (of interest), and $y$ is our target matrix (dependent values).\n",
    "\n",
    "For a Cholesky decomposition to work we assume that $A$ can be broken up into a product of a lower triangular matrix, $L$ and the transpose of the same matrix, $L^{T}$.\n",
    "\n",
    "Note that this is when $A$ is square.  Of course, with an over determined system, $A$ is not square.  So we factor the product $A^{T} \\cdot A$ instead.  We then assume:\n",
    "\n",
    "$$A^{T} \\cdot A = L^{T} \\cdot L$$\n",
    "\n",
    "For more information on the Cholesky decomposition and it's uses, see the following wikipedia link: [The Cholesky Decomposition](https://en.wikipedia.org/wiki/Cholesky_decomposition)\n",
    "\n",
    "Given that $A$ has a unique Cholesky decomposition, we can write our linear regression system as the following:\n",
    "\n",
    "$$ L^{T} \\cdot L \\cdot x = A^{T} \\cdot y $$\n",
    "\n",
    "Then we break apart the system as follows:\n",
    "\n",
    "$$L^{T} \\cdot z = A^{T} \\cdot y$$\n",
    "\n",
    "and\n",
    "\n",
    "$$L \\cdot x = z$$\n",
    "\n",
    "The steps we will take to solve for $x$ are the following\n",
    "\n",
    " 1. Compute the Cholesky decomposition of $A$, where $A^{T} \\cdot A = L^{T} \\cdot L$.\n",
    "\n",
    " 2. Solve ($L^{T} \\cdot z = A^{T} \\cdot y$) for $z$.\n",
    " \n",
    " 3. Finally, solve ($L \\cdot x = z$) for $x$.\n",
    " \n",
    "We start by loading the necessary libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "from tensorflow.python.framework import ops\n",
    "ops.reset_default_graph()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we create a graph session"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sess = tf.Session()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the same method of generating data as in the prior recipe for consistency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create the data\n",
    "x_vals = np.linspace(0, 10, 100)\n",
    "y_vals = x_vals + np.random.normal(0, 1, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We generate the design matrix, $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create design matrix\n",
    "x_vals_column = np.transpose(np.matrix(x_vals))\n",
    "ones_column = np.transpose(np.matrix(np.repeat(1, 100)))\n",
    "A = np.column_stack((x_vals_column, ones_column))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we generate the "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create y matrix\n",
    "y = np.transpose(np.matrix(y_vals))\n",
    "\n",
    "# Create tensors\n",
    "A_tensor = tf.constant(A)\n",
    "y_tensor = tf.constant(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we calculate the square of the matrix $A$ and the Cholesky decomposition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Find Cholesky Decomposition\n",
    "tA_A = tf.matmul(tf.transpose(A_tensor), A_tensor)\n",
    "L = tf.cholesky(tA_A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We solve the first equation. (see step 2 in the intro paragraph above)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Solve L*y=t(A)*b\n",
    "tA_y = tf.matmul(tf.transpose(A_tensor), y)\n",
    "sol1 = tf.matrix_solve(L, tA_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We finally solve for the parameter matrix by solving the second equation (see step 3 in the intro paragraph)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Solve L' * y = sol1\n",
    "sol2 = tf.matrix_solve(tf.transpose(L), sol1)\n",
    "\n",
    "solution_eval = sess.run(sol2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract the coefficients and create the best fit line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "slope: 0.944118894701\n",
      "y_intercept: 0.227194921431\n"
     ]
    }
   ],
   "source": [
    "# Extract coefficients\n",
    "slope = solution_eval[0][0]\n",
    "y_intercept = solution_eval[1][0]\n",
    "\n",
    "print('slope: ' + str(slope))\n",
    "print('y_intercept: ' + str(y_intercept))\n",
    "\n",
    "# Get best fit line\n",
    "best_fit = []\n",
    "for i in x_vals:\n",
    "  best_fit.append(slope*i+y_intercept)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we plot the fit with Matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEACAYAAACnJV25AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX5+PHPCQEhhITFEGRLQtSiUnH5VhFREgS1tqhY\nbdUECSr4UwRkUfRrIcFUv2oDiGBR3BCBtoqtqKgVS4NFBFzLqgIJAUMBAdkEgZDn98fNMpOZSWa5\ns2Tmeb9e8zIz9869Z2bkmTPPPec5RkRQSikVPeLC3QCllFL20sCulFJRRgO7UkpFGQ3sSikVZTSw\nK6VUlNHArpRSUcbrwG6MedEYs8sYs8bhsSeNMRuNMV8ZY94wxiQFp5lKKaW85UuP/WXgqjqPfQCc\nIyLnAZuAh+xqmFJKKf94HdhFZDnwQ53HPhSRyqq7K4HONrZNKaWUH+zMsd8OvGfj8ZRSSvnBlsBu\njHkYOCEiC+w4nlJKKf/FB3oAY8wQ4BqgXwP7aVEapZTyg4gYX/b3tcduqm7WHWOuBh4ArhWRY140\nTm8i5Ofnh70NkXLT90LfC30v6r/5w5fhjguAFcCZxphtxpihwAwgEVhijPnCGPMnv1qhlFLKNl6n\nYkTkVjcPv2xjW5RSStlAZ56GQVZWVribEDH0vail70UtfS8CY/zN4fh8ImMkVOdSSqloYYxBfLx4\nGvComEClp6dTVlYW7mYoP6WlpbF169ZwN0Mp5SDsPfaqb6OQtEHZTz8/pYLLnx675tiVUirKaGBX\nSqkoo4FdKaWijAZ2pZSKMhrY65Genk5CQgLJycm0bduWPn368Nxzz3l1sbCsrIy4uDgqKysb3Fcp\npewU9uGOnpSWljFx4hzKyyvp1CmOwsI8MjLSQnoMYwyLFy8mOzubQ4cOsWzZMkaNGsWqVat46aWX\n6n2uiOiIEaVUeISwkI244+7xkpKtkpk5TuCwgAgclszMcVJSstXtMdyx4xjp6enyz3/+0+mx1atX\nS1xcnKxfv14WL14s559/viQlJUnXrl2loKCgZr+uXbtKXFycJCYmSqtWrWTlypWyZcsW6devn7Rr\n105SUlIkJydHDhw44HV7IpGnz1UpZY+qf2O+xVtfn+DvzZfAnpNT4BCQpSYw5+QUuDmCe3Ycw11g\nF7GC9rPPPivLli2TdevWiYjI2rVrpUOHDrJo0SIREdm6davExcVJZWVlzfM2b94sH374oZw4cUL2\n7Nkjffv2lTFjxnjdnkikgV0pe5SUbJWcnALJypokOTkFNZ1QfwJ7RKZiyssrgZZ1Hm3Jjh3e56vt\nOIYnHTt2ZN++fVx++eU1j/Xo0YObb76ZZcuWce2119Y8LlUpGYDMzEwyMzMBaNeuHWPGjOGRRx4J\nuD1KqcattLSMAQNmsGXLZKy49SMrV+azZMlIv44XkRdPO3WKA36s8+iPdOzofXPtOIYn5eXltG3b\nltWrV9OvXz/at29P69atee6559izZ4/H533//ffccsstdO7cmdatW5Obm1vv/kqp2DBx4hyHoA7Q\nki1bJjNx4hy/jheRgb2wMI/MzHxqA/OPZGbmU1iYF9JjuPPpp5+yY8cO+vTpw6233sr1119PeXk5\n+/fv56677qq5WFrdS3f00EMPERcXx7p169i/fz/z5s3Ti6tKKdszDBGZisnISGPJkpFMnFjEjh2V\ndOwYR2HhSJ9GtNhxDEfVo2Luu+8+Bg8ezDnnnMPhw4dp06YNTZs2ZfXq1SxYsICrrroKgJSUFOLi\n4tiyZQtnnHFGzTFat25NUlIS5eXl/PGPf/SrLUqp6FKbYXAM7gFkGHxNyvt7w4eLp5EiPT1dEhIS\nJCkpSVq3bi29e/eWWbNm1VwQfeONNyQtLU2SkpJk4MCBMnLkSBk8eHDN8/Pz8yUlJUXatGkjq1at\nkvXr18uFF14orVq1kvPPP1+mTp0qXbp0CdfLs0Ukf35KNRb1jeLDj4unWt1RBUQ/P6XsUT3vpjbD\nYM278ae6owZ2FRD9/JQKLi3bq5RSSgO7UkpFotLSMnJzJ/v1XE3FqIDo56eU/ZwnLCVqKkYppRo7\n1wlLvvE6sBtjXjTG7DLGrHF4rI0x5gNjzDfGmH8YY5L9aoVSSqka7icsec+XHvvLwFV1HnsQ+FBE\nfgYsBR7yuyVKKaUATyVRvOdTjt0Ykwa8LSLnVt3/GugrIruMMR2AYhHp7uG5mmOPQvr5qVhlx5oR\n9R07nDn29iKyC0BEdgIpAR5PVZk1axYdOnQgKSmJffv20apVK7Zu3er18+Pi4igpKQHg7rvv5tFH\nHw1SS5UKrerRItnZ+eTmTqa0tCwsbRgwYAbz54+nuHgy8+ePZ8CAGW7b4k97q0ui5OQU+dW+QHvs\n+0SkrcP2vSLSzsNzJT8/v+Z+VlYWWVlZEd3jS09PZ/fu3cTHx9O0aVN69+7Ns88+S6dOnQI6bkZG\nBi+++CL9+vVzu72iooKkpCRWr15Njx49XLYPHTqULl261Fvyt0mTJmzatIlu3boF1NaGRPLnp6KP\nu/K2mZlWeVu7est1z+fYKx8+vD+zZ3/IkiX/YffuV6lb2yUnp4h58/Kdnu9re4uLiykuLq65P3ny\nZJ977L7We0kD1jjc3wikVv3dAdhYz3Pd1kjw9HgkSE9Pl6VLl4qIyLFjx+T222+XQYMG2XJcdwt4\nVNu+fbvExcVJRUWF2+15eXkyceLEes9hjJEtW7YE1E5vRPLnp6KPHQvoeMu1fssGiY8fUnV/Up02\nWLfs7Em2txc/asX4mooxVbdqbwF5VX8PARb5eLyIJ1W90WbNmnHjjTeyYcOGmm3Hjx9n/PjxpKWl\ncdppp3HPPfdw7NgxAPbu3cvAgQNp06YN7dq1o2/fvgDcdtttbNu2jYEDB5KUlERRkfNPrU2bNtG9\nu3WZok2bNvTv3x+oTa08//zzzJ8/nyeffJKkpCSuu+66Bl/D0KFDmTRpEgDLli2jS5cuTJ06ldTU\nVDp16sScOXO8ek1KhVswF9Cpy3XI4WtUVDxTdd+79R5C2V5Hvgx3XACsAM40xmwzxgwFHgcGGGO+\nAfpX3bePMfbeAnDkyBH++te/cskll9Q89sADD7B582bWrFnD5s2bKS8vr0mPTJkyhS5durB37152\n797NY489BsDcuXPp2rUr77zzDgcPHmT8+PFO5znjjDNYv349AAcOHODDDz+seius9g8bNoycnBwe\neOABDh48yKJFvn+X7ty5k0OHDrFjxw5eeOEFRowYwYEDBxp8TUqFWzAX0KnLNSg73s8DGl7vIZTt\ndeJrF9/fG/6kYtz91gnk5qP09HRp1aqVtGnTRuLj46VTp041a5yKiLRs2VJKSkpq7q9YsUIyMjJE\nRGTSpEly/fXXy+bNm90et75UTPV6qSdPnqx5zDG14msqxnH/4uJiSUhIcDp2+/btZdWqVQ2+Jnfq\n/fyUspkdi9R7yzWNUvf+VoHfS2rqYKc1SqvbmZNTIBdfPFoSE4f63t6KCpG//12kT5/oWfM0kixa\ntIjs7GxEhDfffJPLL7+cjRs3YozhyJEjXHjhhTX7VlZW1qRu7r//fgoKCrjyyisxxjBs2DAmTJgQ\nrpfhpF27dsTF1fYYEhISOHz4MN9//329r0mpcLN7AZ36FBbmsXJlvkM65rfEx49wSMecSmbmUZYs\nKXQ6v+sF040kJt5Cjx7dycxsWX97jxyBOXNg2jTYvNnvtkd2YI+AgFId1IwxDBo0iLvuuovly5cz\naNAgEhISWL9+PaeddprL8xITEykqKqKoqIiNGzeSlZXFRRddRHZ2tttl83wR6PM9OfXUU+t9TUpF\ngoyMNKeRJ47sHFvu7ktk+PBhzJ5d/5eKa27+LA4f/jOZmUUe283OnTBzJsyaBfv2+dVeR5Ed2CPM\nokWL2L9/P2effXZNL/y+++5j5syZpKSkUF5ezvr167nyyitZvHgx3bt3JzMzk8TEROLj44mPt97u\n1NRUSkpKPA53BOrtJVc/324NvSalIpm7oYUrVwY2FNLdl8jll19a73N8umC6bh1MnQrz58Px487b\nWreGu+6CJ57wud1aBKwB1aNXkpOTmThxInPnzq0ZtfLEE09w+umn06tXL1q3bs2VV17Jt99+C1ij\nW/r370+rVq249NJLGTFiBJdddhlgLWpdWFhI27ZtmTp1qtvz1u2VO96/4447WL9+PW3btuWGG27w\n6vn1cdz38ccf9/ialIpkrj3llmzZMpmJE+eEtB0NXjAVgSVL4Oqr4ec/h5dfdg7q6ekwfTps3w6P\n+zceRcv2qoDo56ciRXZ2PsXFrvXLs7PzWbp0clBLADjyOClp8V1krFxh9dDXrHF9Yq9eMG4cXH89\nxNcmU/xZQUlTMUqpqFDbU3aeDdqxY1xQ0jSe1M3Nn97uGI9nNKFtdl/473+ddzYGBg2yAnrv3vY1\nwtdhNP7eaIQzT1XD9PNTkaK+oZChnLFaY8sWkZEjRVq2dB16nZAgMmKEyKZNDR4GHe6olIpV9Q2F\nDOkM0E8+sdItf/sbVNY5/mmnwciR1kXRtm3dP98GGtiVUlHD01DI+tI0tjh5EhYtgilTYMUK1+0/\n/zmMHQu33AKnnGLPOeuhF09VQPTzU41B0KpCHj5sjWp56ilwNwT5yiut/PmAAX6XNfHn4qkGdhUQ\n/fxUY1E9KmbHjkqSkg4iEs/Bgwn+jZDZscOaUPTss/DDD87bmjaFnByrh/7znwfc7kYZ2NPT0ykr\nC32hfGWPtLQ0nxYAUSrcAuq9r1lj5c8XLIATJ5y3tWkDd98N995r5dK9bEtDQzD9CexhHxWjlFKh\n5PMImcpKkX/8Q2TAAPfFBbt1E5kxQ+TwYZ/a4W1BM0JQj10ppRo1r0fIHDtmFeTq2ROuusqaLeqo\nd2944w349lurl96y7jHdq14qr1evMUGbKaujYpRSMaXBETL79lm58xkzrOJcjuLi4IYbrAuivXr5\nfG7nNNCTBGsIpvbYlVIxpbAwj8xM10Uy/u+OLKvn3aULPPywc1Bv2RJGjYJNm+D11/0K6lC3nk3w\nFuHQHrtSKqY4TWQqP0nfpuWM+GkzbftNcd25Y0croA8fbl0cDZBzGigPaxUm54u4hYUjAz6PBnal\nVMzJ6NKJededbU0oKl7lsn1DsxRS/u8hUu4dAc2a2XZe5zRQGjASeJzU1DL698+0bdEQTcUopWLH\noUNWSdwzzoDf/hZWOQf197iaK/iQc46XMOaLg7YGdXCXBrJWYfrkk0Lmzcu3rSCZ9tiVUtGvvBye\nfhqeew6qFm6vdtw04VUZwlTGsoFzah4PRh2ZUC3tp4FdKRW9/vMfK93y5z9DRYXztrZt4e67GbPx\nJ/70N8dhh2BrHZk66lvazy5hn3mqlFK2EoH33+foHx6lxYqPXbeffjqMGQNDhkDLlsGrI2OTRllS\nQCmlbPHTT9baoVOnwoYNLps/bd6Jrk9NIvXOO6BJE6dtjnVkrPRIcFZX8kfYArsxZgxwB1AJrAWG\nisjxOvtoYFdK2W/vXpg1yyrKtWuX06aTxLGQG5nCOD7lHHJyioKeBrFbWJbGM8Z0xBqz011Ejhtj\n/grcDMwN9NhKqdjh85qkmzbBtGnWtP+jR502HWnSjOdO3sN0RlNGes3jQVlYIwLZdfG0CdDSGFMJ\nJAA7bDquUioGeL0mqQgsX25dEH3rLeu+o06dYPRo7vt0L8+/PpFQXRCNOL5WDXN3A0YBh4BdwKse\n9vGp8plSKnY0WHHxxAmRv/xF5Be/cF9h8bzzRObNEzl+XES8r5zYGBCONU+NMa2B67CmUR0AFhpj\nbhWRBXX3LSgoqPk7KyuLrKysQE+vlAqQzymQIPBUcfGHbUetdMv06eBu3YZrrrEKcmVnO61QFKrx\n4tXsfA+Li4spLi4OrEG+fhPUvQE3As873B8MzHSzX3C/1pRSPouUnm3dHnsnPpEnuEQOmKauvfNT\nThG5806RDRtqXkNOToFkZU2SnJyCkLc92O8hfvTY7QjsF2GNhGkOGGAOMMLNfra8SKWUfXxedCJI\nqoPjeXwsr3K9HCfOJaBXtG0rMmmSyM6dLs8L5xdTsN9DfwJ7wFcSRGQ1sBD4EvhPVXCfHehxlVLB\n5/WiE8FUWUnGhnWsbb+CL7mUXN6kKbXn/4YzuYvp3Dngbpg8GVJTa7Y5l8G12m7XYhXeioj3sA5b\nRsWIyGSs2pNKqUakwUUngumnn+DVV60c+saNtKizuZi+TGEci/kVwnZSi8eQnZ3vlMOOhKAa1vfQ\nE1+7+P7e0FSMUhEnLKmM3btFCgpEUlJc8+dNmsjHaT3kQj5yeHirwBi3bYyEVFIk5ti1pIBSMcLT\nyI2QTaf/5htruv/cuVZv3VGrVjBsGIweTelJqTOmfSLwIHV7xDk5RRQW5kVEnZdgvodaK0Yp5VbY\nCl2JwEcfWROK3n7bdXuXLjB6NNx5JyQnO7W3OlCuX1/K7t2uE9mzs/NZunRyRNd5sYM/gV1TMUrF\ngJCnLI4fF1mwQOTCC13TLSBywQXW9qoJRZHU9nAPn6yLcExQUkpFvpBdZDxwAF54wZpQtH27y+Yf\nr+jPNNOZf57oQqfF31LYa0eDvevCwjxWrsx3+bVhx9qgdXld2iDCaWBXKgYEfeTGtm1WMH/+eWv5\nOUfNm8OQIXx302/Juutdn4NmKGeReh4+2ciqQvraxff3hqZilAqboI3c+PRTkZtvFmnSxDXdkpJi\njX7ZvVtEImcyVH2ysia5zRxlZ08KW5vQVIxSyh1be72VlbB4sXVBdNky1+3du8PYsZCbCy1qR6dH\nwpjzhkTkmHQ/aGBXKkY4rrXpV9Gqo0etoYrTpllDF+vKzrYC+jXXUFq2nYnDnqS8vJLk5IOIxLNh\nwxYiPWiGMp8fVL528f29oakYpSKCz2mZXbtE8vNFTj3VNUcRHy9y660in3/u4fiOk4s8TzSKJNWj\nYrKzG++oGA3sSsUYr3PdGzZYVRRPOcU1oCclidx/v8i2bQ0cv+65tgr8XlJTB0dE0GwM/AnsmopR\nKsbUm+sWgX/9y8qfv/uu65PT0qwJRXfcAUlJXhy/7rnSgELOPju/cY0yaWQ0sCsVY9xdIIxnPzcd\nWwsXXghffun6pF/8wlrQ4je/gfj6w4bz8aPjYmRjoyUFlIoxjpNwkjnBMJ5hbJPHOe3kYecdjYFr\nr7UCep8+TisUeXt82ANMBwoJZy2XxkxrxSilvLL938tZN/x++m76goSTx503tmgBeXlw332UNj3F\nryXfHOu3JCVZo2IOHUqIylouwaaBXSlVv08/haIiWLjQGo/uqH17uPdeuPtuOPXUkBQOi4T1ViOd\nFgFTSrmqqBD5+99F+vRxHd0CImefLfLCCyJHjzo9LdgzRSNhWbvGgHAsjadUrCotLSM3dzLZ2fnk\n5lrlYyPKkSMwaxacdRYMGgTLlztvv+IKa+TLunXWKJfmzYHa1/XOO1sI5kzRSFjWLlrpqBil/BDR\nVQB37oRnnoE//Qn27XPeFh8PN99sXRA97zyXpzq/riKCOaKlMZQYaKy0x66UHyKyt7l+vdXzTkuD\nP/zBOagnJ8MDD0BpqbXOqJugDnVfVx6QjxXcoXZ6fZ7L8/z59VI7LNKRDoW0ha+5G39vaI5dRZGI\nqQJYWSmyZInI1Ve7z5+np4s89ZTIwYNeHc71dW0VKJDWrT3PFPU3V645du+gJQWUCg3XC4vWVPn2\n7UM0Vf7YMZFXXhHp2dN9QL/oIpHXXhM5ccKnw/pzwTSQi6yRVpclEmlgVypEPBe6sqfn6XF5tn37\nRB5/XKRjR9dgbozIoEEiy5dbPfmAX5d3ryVifr1EqbAFdiAZeB3YCKwHLnazT7Bfv4oRkbImZXU7\nUlMH2Tos0F1w7dv1TtmflyfSsqVrBE1IEBkxQmTTJltfV0O96Or92re39/UrZ+EM7HOAoVV/xwNJ\nbvYJ8stXsSAS87J291gdUxsX84m8xo1SQZzrCVJTRR59VGTvXptfUcOC/YtF1fInsAc83NEY0wq4\nTETyqqJ3BXAw0OMq5U6o16R0nBlZvWDEwYMJTrMk7V5157/fVTCIfzCOKVzKCpftx8/8GXPadmdh\ns3Nov+EEhQcOkdG2rV/n8pfz59ASGA08TmpqGf37ZwZtTVLlJV+/CeregJ7AKuBl4AtgNtDCzX5B\n/l5TsSCU+Vxve6UN/YrwOnV0+LDIjBmyM7GN6wsEWdOhm/x3ziuS2W1s2HvHmlcPHcJUjz0euAAY\nISKfGWOeAh7EGgDrpKCgoObvrKwssrKybDi9iiXBWJPSU70S515pEbUVCqHuLwVP64l6NZHpv/+F\nmTOtWaI//ECqQ9uO05QF/JbXOzdh5kePMH7iHLaUPOKxHaESLWuDRqLi4mKKi4sDO4iv3wR1b0Aq\nUOJwvw/wtpv9gvu1pmKCXTn26l70xRePlsTEoW6P59wr9a+HWu9QwLVrRfLyRJo1czlwRXKyLDqn\nj9zYe4xTL9+XcebBvMgcidc6ohXh6LGLyC5jzHZjzJki8i1wBbAh0OMq5U5GRprH3rG3XKfNP4q7\nHrAdC0a4TpsX+rOCUe+9CvMLXJ/QrRuMGUOTvDyuTUzk2jqbndtUBswAJrN/f0vmz6/9NQDYXvKg\n7i+bl14axOzZ/n8OKoh8/SZwd8PKs38KfAX8DUh2s09wv9aU8pJzL9pzT9yOkR/V52rKMbmNOfIV\n57rNn8sll4gsXGhVYqyHc5s8/xqwuzKj9tDDhzDl2BGR/wC/sONYSgWbcy/ac0+87q+DpCRBpMBh\nwYiGe6iPjrue8967mlv3baYjO522VRrD0auupuWkiXDJJV613bFNixdvYf9+90W0rL6UfQW2Qj0a\nSQVGqzuqmOOczsjDus7vvJhEYaGVzsjISPMvcG3ZAk89RdpLLzH+yBGnTYdpyUvczlMynLhNc1jS\noSMZPhy6uk25uZOZP7++9JB9Fze1EmMj42sX398bmoqJapEyG9QbrmmFDZKYOFB69bo/8LZ//LHI\nDTdY0/vrpFv2tUiUCUyW1uzzOT3i7v2tLz1id+ok2ItuKM/QWjEqHBpj/tXW4lMVFSKvvy7Sq5f7\n/Pm554q88or0v/xhv0bWNBTAPb0OO19jY/yMo4UGdhUWMdubO3RI5OmnRTIy3Af0q64S+eCDmoJc\nDVWEXLZsudtfPb68v8Ee4qiVGENPA7sKi5ibhVheLvLggyKtW7u+6GbNRIYOtcao11H/KJsNEh8/\nxG2P2Nv3V3vV0cmfwK7TxFTAYmYlnDVrYMgQSE+Hxx+H/ftrt7VtCw8/DGVl8NJL0KOHy9OrR7Tk\n5BSRmjoG55msr1FR8Qyuo07meP3+RuSqTiosouxfngqHwsI8MjPz8WYJtUZHBN5/HwYMgJ49Ye5c\nOHGidntmplUOYNs2azm6Dh3qPVz1iJazzvo5zqNMPI868fb91ZErqpoOd1QBs2M2qL881XkJ2LFj\nMH8+TJ1qrSVaV58+1oLQAwdCkyY+H9611or34+k9vb9av0XV8DV34+8NzbErmwUlp7xnj0hhoVXr\nvG5SOy5O5KabRFauDELbPefY/T+m5tijAX7k2I31vOAzxkiozqUah0B729YEnfHU7aHm5PgxG3LT\nJpg2DebMgaNHnbe1bAl33AH33QcZDU8l8vZ1Ve9X3QsfPrw/s2d/6NAr9/3XR91j2vYLRoWNMQYR\nMT49yddvAn9vaI9dObCjdxnwaJzKSpF//1vk+uvdTiiSjh2t9UX37Qvp61LKEToqRjUWdozgcD9a\nZCOlpevIzram3JeWlrk+saICXnsNevWCyy6DN9+0YnC16oukpaUwYQK0aRPS16VUoPTiqQoLO0Zw\nFBbmsXJlvkMg3Uh8/BNs3TqXrVvdlKo9dMgaivjUU7B1q+sBf/lL64Jov35gfPvla+frUipQGthV\nWNgxgqPuaJHS0nVs3TqXur3laeMm8vTp8TB7Nhw44HyQZs1g8GAYOxbOPjuwF4WOTFERwtfcjb83\nNMeuHAQjF103596TL2UuuXLCxLnmz9u1E5k4UWTnTpd2BTIlX3Psym7oqBjVmNg9giM3dzIL5o/l\nav7NOKZwBUtddzrjDKt3ftttkJDg0p66qw5lZvq+6pCOTFF28mdUjAZ2FVJBm1D00098/9R0Dkx6\nktNP7HPdfvnlVv7817+GOPdpEVuHTyplE38Cu+bYVci46xEHug4ne/bArFkwcyYpu3eT4rDppDEc\n/dWvSJw0CX7R8AJfeuFTRQu9oqPqVVpaRm7u5PqHD3rJ1qGA334Ld98NXbvCpEmwe3fttlatrAWh\nS0tJfPttr4I6xFAxMxX1tMeuPLK7hx1wj1gE/v1vmDIF3n7beew5QOfOMHo0DBsGyck+t891+KTz\nMnlKNRYa2JVHdi9g7PdQwIoKWLjQCuiffea6/YILrPz5TTdB06Y+t6taOIuZKWUnDezKI7tzzj73\niA8ehBdfhOnTrTrndf3611ZA79vX7YQify7U+r14tVIRRAO78sjuyTZe94i3b4enn7YmFB086Lyt\neXNrqOKYMdC9u8dzBeVCrVKNha8D3z3dsC7EfgG85WF7cEbvq6AJ+WSbzz8XufVWkfh41wlFKSki\nBQUiu3d7daiYXYdVRR38mKBkZ499NLABSLLxmCqMQpJzrqyEd9+18ufFxa7bf/Yza0LR4MHQooXX\nh9WhiyqW2RLYjTGdgWuAR4GxdhxTRYag5ZyPHoV586wVir7+2nV7VpaVP7/mGo8TiuqjNVtULLPr\n//JpwP2ATi1V9du9GwoKIC0Nhg93DupNmsCtt1ojX/71r3pniTYkqtdhVaoBAffYjTG/AnaJyFfG\nmCzA49TXgoKCmr+zsrLIysoK9PTKQdCm69vQprhvd5P3wyqytq0l7vhx552SkqwgP2oUdOliy3l1\n6KJqrIqLiyl2l5b0QcC1YowxjwG5QAXQAmgF/E1EbquznwR6LuWZXQWsbG1TyVZ+f9kEbt5xiIG8\n57pDly7WcnN33mkFd6WUi7AXATPG9AXGici1brZpYA+iSCpgVfrtZt4dOoE+q5fSs2K/y/aSth3p\n9swU+M37MgXaAAASJUlEQVRvAppQpFQs0CJgMSwiRoEcOMDeJ56k+ZNPM+LkYZfNbzGQIsYTf+6H\nLL355tC1S6kYY2tgF5FlwDI7j6m8E9ZRINu2WbNDn3+edocOOW06SnNeYQjTGMO3/Az4kZxO/wp+\nm5SKYVqPPUrYvUiEVxdgP/vMGn/++utw8qTTpl20Zya5zOIoe/ljQG1SKpaFPcde74k0sAddoCv3\nePXlUFkJ77xjBfSPPnI5RnnSqeQfLGAed3CM5kAZ8AKpqWX075/p10idSBzto1So+BPYdc1TVaPe\nafg//igya5bImWe6TvcHkX79RBYvlpLNJbaWIdA1RFWsI8wlBVQE86bX6+4CbHsOc+Xyf0LXGbB3\nr/NB4+Phd7+zZoiefz4AGWDr+HG7Sgdrr1/FEg3sMcDbSoeOF2DPYgNjmUou82hedsz5gMnJtROK\nOnd2OZ+/ZQjcBV87RvtopUcVc3zt4vt7I8pSMSUlWyUnp0CysiZJTk5BRKcGvK10WLKlVHJPu1EW\nc6X7dEtamsi0aSIHD9reRk8pl2uvHR9wlUat9KgaMzQVExqNrQfYUK+39JtNvD90Av3XfsKrh3e6\nHuCii6x0yw03WOmXIPCUcunR4/dkZga2XF1EjPFXKoQ0sPvB7iXjgs3TGPfTTz3GvgcfokXRTO6u\nM6FIjMFcd50V0C+91O0KRXbyFHwPHkxiyZLbA8rZa6VHFWs0sPuhsfUA6y5Jl8YG8pMHM+Tdb4j7\n8UenfY/QgjnksOXXLZjy96dD1sb6gm+gpYN1kWoVazSw+6Gx9QCrKx3Oufs+sr9YweV7NhJ3wHlO\nwU5SmcFInuX/sY92ZB8O7S+PYAZfrfSoYo1OUML3oXCRWEnRo5Mn4a23rAlFH3/ssnl7cgqTDjzC\nAoZynFOqHg1T8bAAJ1gpFY105qkf/A3SER+EjhyBOXNg2jTYvNl1e//+MG4cpWd2Z8CVMxvHl5RS\nMUgDux8iqdytLXbuhJkzYdYs2LfPeVt8PNxyi3VBtGfPmocj/ktKqRimZXv90NguhHq0bp21fuj8\n+VBnhaJD8c0pPvNCer4wha6XXOzy1KCta6qUCouYD+y+XAiNuGnpIvDPf1r58/ffd9m8LT6ZooqH\neanibn7cYMgcnM+SJR38bnMoX3/EvddKNSa+zmjy90aYZ556minqbZEpu4pR+TNj1eU5X38r8sor\nIuee636G6MUXy1OX3SRxHHA729LfNoSqGJcW/lKqFn7MPI2JwN5QoKgOdNnZrkG/OgCmp98Q8LR0\nfwKW43Nas08mMFn+2yTRNZgbI3LDDSIffywiIllZk9zG/F69RvsVNEM5LV9LAChVSwO7B/4ECtcg\n/LDbQJmdPSmo7cjJKZAM1sp0RsohWro2ICFB5N57RTZv9upc/n5Befqi8OX1eyuU51Iq0vkT2CNz\nRo3N/LlA6lo2oClWLt6Rb5OSfG7HJ58w7B9/ZRM9GcUMEh3Ov6dZIjz2GGzfDjNmQGam01MLC/PI\nzMx3aLM1jLFDh9N9a0OV2msRjoIzKSuU51IqGsXEvxR/AoVrEM4DXANlYWGey3NLS8vIzZ1MdnY+\nubmTKS0t874dJ0/CG29A797Quzd992ykCbVBdy09GMJzjL/hPnjoIUoPHHJ7rurZljk5RWRn55OT\nU8SSJSPJzEzw+b0Az18U7l5/oEJ5LqWikq9dfH9vRHCO3R33qYwNkp5+g0su3ttz1duOw4dFZswQ\n6dbNNQcB8j5XyJW8L3DIu+PZ+F44PtfdtYhgCOW5lIpkaI7dM18Dhb8BsKE8et12lK1cLfLQQyJt\n2rgE82M0lZfIkx68KYmJA6VXr/ud2u7vRUYNmko1Hv4E9pgZx+7rJBx/C0c1lEevaceaNdaEossu\nhRMnnPbeH9eUZypHM5Ox7OQ068HD/cnMdJ4N6+/kKjsmJOk4c6UiV8CB3RjTGZgLdABOAs+LSOjq\nvQaRPwGw3glPIvDBB9aEoiVLXJ67mW5MYyxzKss4wpN1troG7HBVmWxsC40oFXN87eLXvWEF9POq\n/k4EvgG6u9kvyD9YIoO7FM5Z3UbL7ieeFOnRw23+/OuULnI9CySOiqqHvFzKLkwTeXScuVKhQzhS\nMSKyE9hZ9fdhY8xGoBPwdaDHbowcUziHyn4k59Dn3LBjLfETvnfeMS7OWmpu3Dj+30PvUVx8i8PG\nPKwROPXXJg9XnfGoqa+jVJSyNcdujEkHzgNW2XncxiajsoJ5rb+Hv79slc911LIl3H473HcfdOsG\nQKdO/8A5pZIG3EF6+m1kZPSoN2CHo4BXY1toRKlYY1vZXmNMIlAMFIrIIjfbJT+/NgBlZWWRlZUV\n8Hkj5iKeCKxYYeXP33zTuu+oY0cYORLuugvatHHa1KgW7qDxtVepxqS4uJji4uKa+5MnT0bCUY/d\nGBMPvAO8JyLTPewjdn2JVIuIAFNRwa7nZnMo/zFO31vuur1nT6v++e9+B82aeTxMY6uJ3tjaq1Rj\nFbaFNowxc4E9IjK2nn1sD+xhXSTj8GF46SVOFBXRdPt2l81H+vYl4fe/hyuuAOPTZ6KUUjX8CewB\nJ0WNMZcCOUA/Y8yXxpgvjDFXB3pcb4TlIl55OUyYAJ07w+jRTkH9GM14kds5h9UM75xtLT+nQV0p\nFWJ2jIr5GGhiQ1t8FtKLeF99ZU0o+vOfoaLCadNe2vIn7uEZRrCLDgCk7njH/jYopZQXGt3MU8eL\npcnJB+na9X/Ztu0x6hsW6DcRa2WiKVOslYrqOuMMXm59JiM+fZmjpDhs0BEiSqnwaVSLWbu7WNql\nyxjOPz+ZQ4cS7LuI99NP1tqhU6fChg2u2/v0sS6IDhxI6bbvwn8BVykVtcJ28dSrE9kQ2IN+sXTv\nXpg1C2bOhF27nLfFxcGNN1oB/aKLnDbpCBGlVLD4E9gbVSrGl4ulPo1v37QJpk2DOXPg6FHnbYmJ\ncOedMHo0pKe7fbq/k4QiZgy+UiqqNKrA7u3FUq+KVInA8uVW/vytt1wnFHXqZAXzYcOgdWvbX4sW\n0lJKBY2vxWX8vWFDETBvi17VW6TqxAmRv/xF5H/+x21BLjn/fJF580SOHw+4vfXRQlpKKW8Q7fXY\nvS165S5lk0gll6xaAaefDmVlrge/5hoYN47S9G5MnPQK5S/8wZb0iKd0ixbSUkoFS6MK7OBdPtsx\nZdOJ7xjF0wxnNq03H3De8ZRTYPBgGDsWzjrL9vRIfcfTQlpKqaDxtYvv740Q1mMvKdkqAzvnyqv8\nTo4T75puOfVUkUmTRHbudHqe3emR+o7nby316mXtsrJ0WTulYgHRnoppUGUlvPceGVOm8NZ3/3Ld\nfuaZVu/8ttugRQuXzXanR+o7nj+11PWCq1LKG9ER2H/6CV591RqyuHGj6/a+fa3x57/6lTUe3QO7\n0yMNHc/XYZITJ85xCOoALdmyZTITJ4ag6JlSqtFo3And77+HRx6BtDQYPtw5qDdpArfcAp99BsXF\nMHBgvUEdoLAwj8zMfKxgDLUlCvL8ap7dx9MLrkopbzTOHvs331i981desXrrjlq1soL8qFHQtatP\nh7V7qTm7j6cXXJVS3mg8JQVE4KOPoKgI3nFTObFLF2tC0Z13QnKy/+eJYBGxsIhSKqSis1bMiROw\ncKE1Q/Tzz123X3ihlT+/8UZo2jTwhkY4rUujVGyJ+MBeUrLV+9ooBw7ACy/A9OngZoUifv1rGD8e\nLr9cF7NQSkWtiA/smZnjGk4jbNtmBfPnn4dDh5wP0ry5NVRxzBjo3j0k7VZKqXCK+MAOh/FYcvfz\nz610y2uvwcmTzk9OSYERI+Cee6y/lVIqRjSCsr3OQ/UMLUhbsxGysmDZMtfdu3e3eueDB7udUKSU\nUspVWHrszTnKbcxlDFPozibXnbOzrQuiv/xlg2PPlVIqmkV8jz2F+xlBMvfwAinscdp20hiOXnsd\niZMmwgUXBLUdusCFUiqahTSwb+NFmnPc6bEDJDGb4Twtd3LKuudZ0qYdGUFsg9ZbUUpFu5DmORyD\n+p6EZMbwOF3YzgP8ke/4WVXdkzlBbYPneivBPa9SSoWKLYHdGHO1MeZrY8y3xpgJ9e37n1M6sPvp\nGfzuf0byFBM4RFLVljKgiMWLt5CbO5nSUjeLYdhA660opaJdwIHdGBMHzASuAs4BbjHGuB1kXjhg\nKEkbPqH9yHs5rUs8tcWxyoAZwHj275/L/PnjGTBgRlCCe229FUdab0UpFT0CHhVjjOkF5IvIL6vu\nP4hVGP6JOvs5lRRwznUXAeOp7UmXAS/Qvn0ZAwZk2npxU+utKKUak7BMUDLG/Aa4SkSGV93PBS4S\nkVF19nOpFVM9OmXx4i3s3z+36tHq3rvvgdfb0S5ab0Up1ViEa7ijuxO6/bYoKCio+TsrK4usrCzm\nzcsnN3cy8+dXl6OdQ21QB28Xk/BltIuvC1wopVSoFBcXU1xcHNAx7ErFFIjI1VX3vUrFOHIOyk9i\nBXZn2dn5LF3q+ng168vBMZ0DTiULlFKqEfKnx27HFcNPgdONMWnGmGbAzcBbvhygekGKnJwiUlPX\n4s/FTR3topRSloADu4icBO4FPgDWA38RETcLj9avOj3yySfT/FpOTke7KKWUJSIX2vDn4qaOdlFK\nRaOIL9vr77l0tItSKlZFZWDXnrhSKpaF6+JpUGltF6WU8k3EB3Yd7aKUUr6J+MCuo12UUso3ER8d\nCwvz/Br+qJRSsSriL56CjnZRSsWuqBwVo5RSsSwqR8UopZTyjQZ2pZSKMhrYlVIqymhgV0qpKKOB\nXSmloowGdqWUijIa2JVSKspoYFdKqSijgV0ppaKMBnallIoyGtiVUirKaGBXSqkoo4FdKaWijAZ2\npZSKMgEFdmPMk8aYjcaYr4wxbxhjkuxqmFJKKf8E2mP/ADhHRM4DNgEPBd6k6FdcXBzuJkQMfS9q\n6XtRS9+LwAQU2EXkQxGpXlV6JdA58CZFP/2ftpa+F7X0vail70Vg7Myx3w68Z+PxlFJK+SG+oR2M\nMUuAVMeHAAEeFpG3q/Z5GDghIguC0kqllFJeC3jNU2PMEGA40E9EjtWzny54qpRSfvB1zdMGe+z1\nMcZcDTwAXF5fUPenYUoppfwTUI/dGLMJaAbsrXpopYjcY0fDlFJK+SfgVIxSSqnIEvSZp8aYq40x\nXxtjvjXGTAj2+SKVMaazMWapMWaDMWatMWZUuNsUbsaYOGPMF8aYt8LdlnAyxiQbY16vmuy33hhz\ncbjbFC7GmDHGmHXGmDXGmPnGmGbhblMoGWNeNMbsMsascXisjTHmA2PMN8aYfxhjkhs6TlADuzEm\nDpgJXAWcA9xijOkezHNGsApgrIicDVwCjIjh96LaaGBDuBsRAaYD74rIWUBPYGOY2xMWxpiOwEjg\nAhE5F+sa4M3hbVXIvYwVLx09CHwoIj8DluLFRNBg99gvAjaJSJmInAD+AlwX5HNGJBHZKSJfVf19\nGOsfb6fwtip8jDGdgWuAF8LdlnAyxrQCLhORlwFEpEJEDoa5WeHUBGhpjIkHEoAdYW5PSInIcuCH\nOg9fB7xS9fcrwPUNHSfYgb0TsN3h/nfEcDCrZoxJB84DVoW3JWE1Dbgfa05ELOsG7DHGvFyVlppt\njGkR7kaFg4jsAKYA24ByYL+IfBjeVkWE9iKyC6wOIpDS0BOCHdjdDXGM6X/IxphEYCEwuqrnHnOM\nMb8CdlX9gjG4//8kVsQDFwDPiMgFwBGsn94xxxjTGqt3mgZ0BBKNMbeGt1WNU7AD+3dAV4f7nYmx\nn1aOqn5eLgReFZFF4W5PGF0KXGuMKQH+DGQbY+aGuU3h8h2wXUQ+q7q/ECvQx6L+QImI7BORk8Df\ngN5hblMk2GWMSQUwxnQAdjf0hGAH9k+B040xaVVXt28GYnkExEvABhGZHu6GhJOI/K+IdBWRblj/\nTywVkdvC3a5wqPqJvd0Yc2bVQ1cQuxeUtwG9jDHNjTEG672IxQvJdX/FvgXkVf09BGiwUxjQzNOG\niMhJY8y9WOV944AXRSQWPyiMMZcCOcBaY8yXWCmp/xWR98PbMhUBRgHzjTFNgRJgaJjbExYistoY\nsxD4EjhR9d/Z4W1VaBljFgBZQDtjzDYgH3gceN0YczvWl99NDR5HJygppVR00aXxlFIqymhgV0qp\nKKOBXSmloowGdqWUijIa2JVSKspoYFdKqSijgV0ppaKMBnallIoy/x9ihTW1hgupHAAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f79b42a82b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the results\n",
    "plt.plot(x_vals, y_vals, 'o', label='Data')\n",
    "plt.plot(x_vals, best_fit, 'r-', label='Best fit line', linewidth=3)\n",
    "plt.legend(loc='upper left')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
